/*
 * generateDag.cpp
 *
 *  Created on: 13/02/2014
 *      Author: mattia
 */

#include "signatures.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <queue>
#include <math.h>
#include <limits.h>
#include <list>
#include <set>

#define INF INT_MAX
#define ATTEMPTS_LIMIT 10000
#define MAX_CI_DIFFERENCE_ALLOWED 1
#define MAX_CNC_DIFFERENCE_ALLOWED 1

using namespace std;

int selectFromDiscreteUniformDestribution(int from, int to);
void fillBasedOnDUD(int n, int min, int max, int * a);
int measureCI(int n, int * numP, int ** P);

int generateDag(int n, float expected_cnc, float expected_ci, int * numP, int ** P, float & cnc, float & ci){

	float current_cnc, current_ci, it_expected_cnc;
	int attempt = 0, nodo, index, i, j, numAristas, sumMi_remaining, M_nodo, m_nodo;
	list<int> remainingNodes;

	it_expected_cnc = expected_cnc;

	float factor = 1;

	do{

		it_expected_cnc = (current_ci < expected_ci) ? pow(1.05, attempt++) * it_expected_cnc
                                                     : pow(0.95, attempt++) * it_expected_cnc;

		sumMi_remaining = 2 * n * it_expected_cnc;

		for (i = 0; i < n; i++){
			remainingNodes.push_back(i);
		}

		numAristas = 0;

		for (i = 0; i < n; i++){

			//Elegir nodo
			index = selectFromDiscreteUniformDestribution(0, remainingNodes.size() - 1);

			list<int>::iterator it = remainingNodes.begin();
			for (j = 0; j < index; j ++){
				it++;
			}

			nodo = *it;

			remainingNodes.erase(it);

			//Elegir numero de predecesores
			M_nodo = (nodo > n / 2) ? (n / 4) : nodo;
			//M_nodo = nodo;
			M_nodo = (M_nodo < sumMi_remaining) ? M_nodo : sumMi_remaining;

			sumMi_remaining -= M_nodo;

			numP[nodo] = selectFromDiscreteUniformDestribution(0, M_nodo);
			numAristas += numP[nodo];

			if (attempt > 1)
				delete [] P[nodo];

			P[nodo] = new int[numP[nodo]];

			//Crear las adyacencias
			fillBasedOnDUD(numP[nodo], 0,  nodo - 1, P[nodo]);

		}

		current_ci = measureCI(n, numP, P);
		current_cnc = numAristas / n;

	}while((abs(expected_ci - current_ci) > (expected_ci * MAX_CI_DIFFERENCE_ALLOWED))
		                    && (attempt <= ATTEMPTS_LIMIT));


	if (false && abs(expected_ci - current_ci) > (expected_ci * MAX_CI_DIFFERENCE_ALLOWED)){
		return -1;
	}else{
		cnc = current_cnc;
		ci = current_ci;
		return 0;
	}

}

int generateDag2(int n, float expected_cnc, float expected_ci, int * numP, int ** P, float & cnc, float & ci){

	float current_cnc, current_ci;
	int attempt = 0, i, j, index, numAristas, availablePredsCount;

	float factor = 1;

	bool haveSucc[n];
	memset(haveSucc, 0, sizeof(bool) * n);
	int haveSuccCounter = 0;

	set<int> availablePreds;
	set<int>::iterator it;

	P[0] = NULL;
	numP[0] = 0;

	do{
		attempt++;

		//Inicializar conjunto de predecesores disponibles
		availablePreds.clear();
		for (i = 0; i < n - 1; i++){
			availablePreds.insert(i);
		}

		numAristas = 0;

		for (i = 1; i < n - 1; i++){

			for (it = availablePreds.begin(); it != availablePreds.end(); it++){
				*it;
			}

			//Contar los predecesores disponibles para i
			availablePredsCount = 0;
			for (it = availablePreds.lower_bound(i); it != availablePreds.begin(); it--){
				availablePredsCount++;
			}

			//Elegir numero de predecesores
			numP[i] = (i < availablePredsCount) ?
					selectFromDiscreteUniformDestribution(0, round(factor * i)) :
					selectFromDiscreteUniformDestribution(0, round(factor * availablePredsCount));

			if (attempt > 1)
				delete [] P[i];

			if(numP[i] > 0){

				P[i] = new int[numP[i]];

				//Obtener indices del set que representan los predecesores
				fillBasedOnDUD(numP[i], 0,  availablePredsCount - 1, P[i]);

				//Crear las adyacencias y remover predecesores seleccionados de los disponibles
				it = availablePreds.begin();
				index = 0;

				for (j = 0; j < numP[i]; j++){

					while(index < P[i][j]){
						it++;
						index++;
					}

					P[i][j] = *it;

				}


			}else{

				numP[i] = 1;
				P[i] = new int[1];
				P[i][0] = 0;

			}

			for (j = 0; j < numP[i]; j++){

				availablePreds.erase(P[i][j]);

				if(!haveSucc[P[i][j]]){
					haveSucc[P[i][j]] = true;
					haveSuccCounter++;
				}

			}

			numAristas += numP[i];

		}

		//Definir predecesores de la actividad sink: n-1
		if (attempt > 1)
			delete [] P[n - 1];

		numP[n - 1] = n - 1 - haveSuccCounter;
		P[n - 1] = new int[numP[n - 1]];
		j = 0;

		for(i = 0; i < n - 1; i++){
			if(!haveSucc[i]){
				P[n - 1][j] = i;
				j++;
			}
		}

		//Ajustar factor
		factor = (expected_ci < current_ci) ? 0.95 * factor
				                            : 1.05 * factor;

		if ((factor > 1) || isnan(factor)){
			factor = 1;
		}

		//Medir complejidad del dag producido
		current_ci = measureCI(n, numP, P);
		current_cnc = numAristas / n;

	}while(((abs(expected_ci - current_ci) > (expected_ci * MAX_CI_DIFFERENCE_ALLOWED))
		   && (abs(expected_cnc - current_cnc) > (expected_cnc * MAX_CNC_DIFFERENCE_ALLOWED)))
		   /*&& (attempt <= ATTEMPTS_LIMIT)*/);


	if (false && abs(expected_ci - current_ci) > (expected_ci * MAX_CI_DIFFERENCE_ALLOWED)){
		return -1;
	}else{
		cnc = current_cnc;
		ci = current_ci;
		return 0;
	}

}

int generateDag3(int n, float expected_cnc, float expected_ci, int * numP, int ** P, float & cnc, float & ci){

	float current_cnc, current_ci;
	int attempt = 0, i, j, k, l, index, numAristas, availablePredsCount;

	float factor = 1;

	bool haveSucc[n];
	memset(haveSucc, 0, sizeof(bool) * n);
	int haveSuccCounter = 0;

	int * P_aux;
	int P_aux_len, pred1, pred2;

	P[0] = NULL;
	numP[0] = 0;

	do{
		attempt++;

		numAristas = 0;

		for (i = 1; i < n - 1; i++){

			//Elegir numero de predecesores
			numP[i] = selectFromDiscreteUniformDestribution(0, round(factor * i));

			if (attempt > 1)
				delete [] P[i];

			if(numP[i] > 0){

				P_aux_len = numP[i];
				P_aux = new int[P_aux_len];

				//Obtener P_aux_len predecesores en el intervalo [0..(i-1)]
				fillBasedOnDUD(P_aux_len, 0,  i - 1, P_aux);

				//Eliminar predecesores redundantes
				for(j = 0; j < P_aux_len; j++){

					pred1 = P_aux[j];

					for (k = j + 1; k < P_aux_len; k++){

						pred2 = P_aux[k];

						for (l = 0; l < numP[pred2]; l++){
							if((pred1 == P[pred2][l]) && (P_aux[j] != -1)){
								numP[i]--;
								P_aux[j] = -1;
								break;
							}
						}

					}

				}

				P[i] = new int[numP[i]];
				j = 0;

				for (k = 0; k < P_aux_len; k++){

					if(P_aux[k] != -1){
						P[i][j] = P_aux[k];
						j++;
					}

				}

				delete [] P_aux;


			}else{

				numP[i] = 1;
				P[i] = new int[1];
				P[i][0] = 0;

			}

			for (j = 0; j < numP[i]; j++){

				if(!haveSucc[P[i][j]]){
					haveSucc[P[i][j]] = true;
					haveSuccCounter++;
				}

			}

			numAristas += numP[i];

		}

		//Definir predecesores de la actividad sink: n-1
		if (attempt > 1)
			delete [] P[n - 1];

		numP[n - 1] = n - 1 - haveSuccCounter;
		P[n - 1] = new int[numP[n - 1]];
		j = 0;

		for(i = 0; i < n - 1; i++){
			if(!haveSucc[i]){
				P[n - 1][j] = i;
				j++;
			}
		}

		//Ajustar factor
		factor = (expected_ci < current_ci) ? 0.95 * factor
				                            : 1.05 * factor;

		if ((factor > 1) || isnan(factor)){
			factor = 1;
		}

		//Medir complejidad del dag producido
		current_ci = measureCI(n, numP, P);
		current_cnc = numAristas / n;

	}while(((abs(expected_ci - current_ci) > (expected_ci * MAX_CI_DIFFERENCE_ALLOWED))
		   && (abs(expected_cnc - current_cnc) > (expected_cnc * MAX_CNC_DIFFERENCE_ALLOWED)))
		   /*&& (attempt <= ATTEMPTS_LIMIT)*/);


	if (false && abs(expected_ci - current_ci) > (expected_ci * MAX_CI_DIFFERENCE_ALLOWED)){
		return -1;
	}else{
		cnc = current_cnc;
		ci = current_ci;
		return 0;
	}

}

int selectFromDiscreteUniformDestribution(int from, int to){
	return from + rand() % ((to - from) + 1);
}

// Precondicion: n <= max - min + 1
void fillBasedOnDUD(int n, int min, int max, int * a){

	//Creating seed
	int rangeLen = max - min + 1;
	int rangePermutation[rangeLen];

	for (int i = 0; i < rangeLen; i++){
		rangePermutation[i] = i;
	}


	//My Fisher–Yates shuffle adaptation
	int swapAux;
	int index;

    for (int i = rangeLen - 1; i > 0; i--){

    	index = selectFromDiscreteUniformDestribution(0,i);

		swapAux = rangePermutation[index];
		rangePermutation[index] = rangePermutation[i];
		rangePermutation[i] = swapAux;

    }

    for (int i = 0; i < rangeLen - 1; i++){

    	index = selectFromDiscreteUniformDestribution(i,rangeLen - 1);

		swapAux = rangePermutation[index];
		rangePermutation[index] = rangePermutation[i];
		rangePermutation[i] = swapAux;

    }


    //Select the first n elems in asc order
    for (int i = 0; i < n; i++){

    	a[i] = rangePermutation[i];

    	// Insertar burbuja en a de modo que quede ordenado
    	for (int j = i; (j > 0) && (a[j - 1] > a[j]); j--){

			swapAux = a[j - 1];
			a[j - 1] = a[j];
			a[j] = swapAux;

    	}

    }

}


void boolArrayIntersection(int length, bool * set1, bool * set2);
int hopcroft_karp(int n, int * numAdjs, int ** G);

int measureCI(int n, int * numP, int ** P){

	int i,j,k,predIndex,pred, dom, rdom;

	//Step 1: Calcular conjunto de dominantes para cada nodo
	//dom(i) = i union {interseccion( dom(j) ) / j es predecesor de i }
	bool domSet[n][n], rdomSet[n][n];

	memset(domSet, 0, sizeof(bool) * n * n);

	for (i = 0; i < n; i++){

		memset(domSet[i], 1, sizeof(bool) * (i + 1));

		for (predIndex = 0; predIndex < numP[i]; predIndex++){

			boolArrayIntersection( i, domSet[ P[i][predIndex] ], domSet[i] );

		}

	}

	memset(rdomSet, 1, sizeof(bool) * n * n);

	for (i = n - 1; i > 0; i--){

		memset(rdomSet[i], 0, sizeof(bool) * i);

		for (predIndex = 0; predIndex < numP[i]; predIndex++){

			pred = P[i][predIndex];

			boolArrayIntersection( n - pred - 1, rdomSet[i] + pred + 1, rdomSet[pred] + pred + 1 );

		}

	}


	//Step 2: Calcular la matriz de adyacencia transitiva de P
	//tam: transitive adjacency matrix
	bool tam[n][n];
	memset(tam, 0, sizeof(bool) * n * n);

	for (i = n - 1; i > 0; i--){

		for (predIndex = 0; predIndex < numP[i]; predIndex++){

			pred = P[i][predIndex];

			tam[pred][i] = 1;

			for (k = i + 1; k < n; k++){

				tam[pred][k] = tam[pred][k] || tam[i][k];

			}

		}

	}


	//Step 3: Calcular el grafo de complejidad
	for (i = 0; i < n; i ++){

		dom = i;

		while ((dom >= 0) && !domSet[i][--dom]);

		for(k = 0; k <= dom; k++){
			tam[k][i] = 0;
		}

		rdom = i;

		while ((rdom < n) && !rdomSet[i][++rdom]);

		for(k = rdom; k < n; k++){
			tam[i][k] = 0;
		}

	}


	//Step 4: Calcular el grafo bipartito equivalente
	list<int> adjacency[n];
	list<int> senderNodes;
	int receiverOrd[n];
	int numReceiverNodes = 0;

	for (i = 0; i < n; i++){

		// Chequear si salen aristas de i y en dicho caso guardar los adyacentes.
		bool sender = false;

		for (j = i + 1; j < n; j++){

			if(tam[i][j]){
				sender = true;
				adjacency[i].push_back(j);
			}

		}

		if (sender){
			senderNodes.push_back(i);
		}


		// Chequear si el nodo i recibe aristas
		j = 0;

		while ((j < i) && !tam[j][i]){
			j++;
		}

		if (j < i) {
			receiverOrd[i] = numReceiverNodes++;
		}

	}


	int size = senderNodes.size() + numReceiverNodes;
	int * bipartiteGraph[size];
	int successors[size];

	i = 0;

	for (list<int>::iterator sender = senderNodes.begin(); sender != senderNodes.end(); sender++){

		successors[i] = adjacency[*sender].size();
		bipartiteGraph[i] = new int[successors[i]];

		j = 0;

		for (list<int>::iterator adj = adjacency[*sender].begin(); adj != adjacency[*sender].end(); adj++){

			bipartiteGraph[i][j] = senderNodes.size() + receiverOrd[*adj];
			j++;

		}

		i++;
	}

	for (j = i; j < size; j++){
		bipartiteGraph[j] = NULL;
		successors[j] = 0;
	}


	//Step 5: Devolver el maximo matching del grafo bipartito
	int ci = hopcroft_karp(size, successors, bipartiteGraph);

	//Liberar memoria
	for (j = 0; j < size; j++){
		if(successors[j] > 0){
			delete [] bipartiteGraph[j];
		}
	}

	return ci;
}


void boolArrayIntersection(int length, bool * set1, bool * set2){

	for (int i = 0; i < length; i++){

		set2[i] = set1[i] && set2[i];

	}

}


bool bfs(int n, int * numAdjs, int ** G, int * match, int * dist);
bool dfs(int u, int n, int * numAdjs, int ** G, int * match, int * dist);

int hopcroft_karp(int n, int * numAdjs, int ** G) {
    int matching = 0, i;
    int match[n], dist[n + 1];

    for (i = 0; i < n; i++){
    	match[i] = n;
    }

    while(bfs(n, numAdjs, G, match, dist)){

        for(i=0; i < n; i++) {

        	if((match[i] == n) && dfs(i, n, numAdjs, G, match, dist)){
                matching++;
        	}
        }

    }

    return matching;
}


bool bfs(int n, int * numAdjs, int ** G, int * match, int * dist) {

	int i, u, v;
    queue< int > Q;

    for(i = 0; i < n; i++) {

    	if(match[i] == n) {
    		dist[i] = 0;
            Q.push(i);
    	}else{
    		dist[i] = INF;
        }
    }

    dist[n] = INF;

    while(!Q.empty()) {

    	u = Q.front();
    	Q.pop();

    	if(u != n) {

            for(i = 0; i < numAdjs[u]; i++) {

            	v = G[u][i];

            	if(dist[match[v]] == INF) {

                    dist[match[v]] = dist[u] + 1;
                    Q.push(match[v]);

                }
            }

        }
    }

    return (dist[n]!=INF);
}


bool dfs(int u, int n, int * numAdjs, int ** G, int * match, int * dist) {

	int i, v;

    if(u != n) {

        for(i = 0; i < numAdjs[u]; i++) {

        	v = G[u][i];

        	if(dist[match[v]] == dist[u]+1) {

        		if(dfs(match[v], n, numAdjs, G, match, dist)) {
                    match[v] = u;
                    match[u] = v;
                    return true;
                }
            }
        }

        dist[u] = INF;

        return false;

    }

    return true;
}
